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Computational Model 

Tracking Primary Electrons, Secondary Electrons, and Ions 
in the Discharge Chamber of an Ion Engine 


Sudhakar Mahalingam and James A. Menart 
Wright State University 
Dayton, Ohio 45435 


Computational modeling of the plasma located in the discharge chamber of an ion engine is an important activity 
so that the development and design of the next generation of ion engines may be enhanced. In this work a 
computational tool called XOOPIC is used to model the primary electrons, secondary electrons, and ions inside the 
discharge chamber. The details of this computational tool are discussed in this paper. Preliminary results from 
XOOPIC are presented. The results presented include particle number density distributions for the primary 
electrons, the secondary electrons, and the ions. In addition the total number of a particular particle in the discharge 
chamber as a function of time, electric potential maps and magnetic field maps are presented. A primary electron 
number density plot from PRIMA is given in this paper so that the results of XOOPIC can be compared to it. 
PRIMA is a computer code that the present investigators have used in much of their previous work that provides 
results that compare well to experimental results. PRIMA only models the primary electrons in the discharge 
chamber. Modeling ions and secondary electrons, as well as the primary electrons, will greatly increase our ability to 
predict different characteristics of the plasma discharge used in an ion engine. 


Nomenclature 


A g 

Grid surface area 

A e 

Magnetic vector potential 

B 

Magnetic flux density 

z 

Axial position or direction 

c 

Speed of light 

A r 

Grid spacing in radial direction 

5 

Electric displacement 

At 

Time step size 

E 

Electric field 

A z 

Grid spacing in axial direction 

H 

Absolute electron charge 

e 

Kinetic energy of a particle 

fB 

Beam ion fraction 

£0 

Permittivity of free space 

H 

Magnetic field strength 

Y 

Relativistic factor 

I 

Current 

X D 

Debye length 

J 

Current density 

P 

Permeability 

k B 

Boltzmann’s constant 

o 

Collision frequency 

1 ,X 

Position 

TO p e 

Plasma electron frequency 

m 

Mass of a particle 

<l> 

Electric potential 

til 

Neutral mass flow rate 


Grid transparency to neutral particles 

N 

Number of particles 

P 

Charge density 

n 

Unit normal 

CT 

Collision cross-section area 

n 

Particle number density 



tip 

Plasma source rate 

Subscripts 

Pnull 

Probability of null collision 

c 

Referring to a constant 

q 

Charge of a particle 

coll 

Referring to colliding particles 

R 

A random number 

e 

Referring to an electron 
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r 

Radial position or direction 

i 

Referring to an i th type charge particle 

S 

Surface area 

ion 

Referring to an ion 

T 

Temperature of a particle 

inc 

Referring to an incident particle 

t,s 

Rotational vectors 

max 

Maximum value 

t 

Time 

n 

Referring to a neutral particle 

u 

Particle relativistic velocity 

pe 

Referring to a primary electron 

_ tip 

Volume of a cathode tip 

r 

Referring to radial position or direction 

V d 

Discharge voltage 

T 

Total value 

V 

Particle speed 

weight 

Particle weighting 

V 

Velocity 

z 

Referring to axial position or direction 



0 

Referring to azimuthal position or direction 


I. Introduction 

Ion engines are one of the most promising spacecraft propulsion devices for future missions to distant planets 
and other celestial bodies. Ion engines have proved themselves to be a successful propulsion system on Deep Space 
1. To enhance the development and design of future ion engines it is necessary to improve upon the computational 
tools available for modeling the performance of an ion engine. The work being presented here is a step in this 
direction. In this paper a computer code called XOOPIC (ref. 1 ) is utilized to model the details of the plasma located 
in the discharge chamber of an ion engine. XOOPIC is a particle-in-cell code developed by the Plasma Theory and 
Simulation Group at UC, Berkeley. While XOOPIC was not developed for use in ion engine modeling it has been 
used for a number of plasma application (refs. 2 through 6). In this work XOOPIC will be used to model the plasma 
in the discharge chamber of an ion engine. 

To see where the present work fits into the current state of discharge chamber modeling a brief history of the 
computational tools that were used for modeling the ion engine’s discharge chamber is given here. In the early 
1990’s Arakawa and Ishihara (ref. 7) developed a complete computational model to predict the behavior of an ion 
engine discharge chamber. Their computational model is inclusive of the models developed by Arakawa and 
Yamada, (ref. 8) Arakawa and Wilbur, (ref. 9) and Brophy and Wilbur (ref. 10). Brophy and Wilbur developed a 
theoretical model to study the performance of an ion engine discharge chamber. Arakawa and Yamada developed a 
computational tool called PRIMA (refs. 8 and 11) to model the primary electron trajectories using Monte Carlo and 
finite element methods. PRIMA tracks primary electrons throughout the discharge chamber and includes the effects 
of magnetic fields; however it does not include the effects of electric fields. Arakawa and Wilbur developed tools to 
model the magnetic field and plasma present inside the ion engine discharge chamber using a diffusion model. This 
model looks at the ions as a continuum and considers the movement of ions as a diffusion process. It does not appear 
as if secondary electrons are modeled. The effects of magnetic fields on the ions are accounted for by using two 
diffusion coefficients, one parallel to the magnetic field and one perpendicular to the magnetic field. Arakawa and 
Ishihara combined these two models to simulate the complete performance of an ion engine discharge chamber. In 
addition they also studied ion motions through the grid holes in their model. Mahalingam and Menart (ref. 11) 
improved PRIMA and upgraded it to handle any shape of axi-symmetric discharge chamber, increased the 
computational speed so much finer discretizations could be used, and altered the particle tracking procedure to some 
extent. Though these models are impressive, there is still a great deal of room for improvement. These models make 
many assumptions like ignoring electric fields and taking all of the primary electrons to be at the same energy. This 
is done to make the computer models manageable. The actual physical phenomena occurring inside a discharge 
chamber is extremely complicated and the need for assumptions will continue to exist. The goal of this work is to 
develop a model that is more realistic than those discussed above. The way to do this is to use a particle-in-cell 
technique for the primary electrons, secondary electrons and ions where all these particles are allowed to interact 
and electric fields are calculated based on the charged particle locations. This is what is done in this work. It should 
be noted that Hirakawa and Arakawa (ref. 12) developed a computational model that includes electric field effects. 
However, their model lacks particle collisions, assumes constant particle energy and the study is limited to the cusp 
region. Recently Wirz and Katz (ref. 13) at JPL developed a 2-D computational tool to model the nonlinear behavior 
of plasma inside the discharge chamber. Their model includes electric field effects, neutral particles, primary 
electrons, secondary electrons and positive ions. Currently their computational model seems to be the most advanced 
computational tool available to model the ion engine discharge chamber. Stueber (ref. 14) at NASA Glenn is 
developing a 3-D computational model to track the primary electron motions inside the discharge chamber. It is the 
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first attempt of a computational tool to include three dimensional effects in the modeling work of an ion engine 
discharge chamber. Stueber’s model only considers primary electrons. The work presented in this paper is the one 
that includes the physical phenomena: tracking of primary electrons, tracking of secondary electrons, tracking of 
ions, interaction of charged particles with neutral particles, electric field effects on particles and the particle’s affect 
on the electric fields. In this work the neutral particles are viewed as a uniform background gas through which the 
charged particles move. 

This paper first presents the models used to determine the static fields in the discharge chamber, both magnetic 
and electric. The magnetic static field is produced by the permanent magnets and the static electric fields are 
produced by potentials applied to walls of the discharge chamber and the cathode. While the static electric fields are 
generated within XOOPIC the static magnetic field results are generated external to XOOPIC. Static magnetic field 
results are passed as an input to XOOPIC. After the static field modeling is presented the mathematical model and 
computational techniques used by XOOPIC to compute what we will call the dynamic field effects and to compute 
the trajectories of the particles are presented. The dynamic electric field as well as the dynamic magnetic field 
modeling is presented in this section because these quantities depend on the location and velocity of the charged 
particles; unlike the static electric and magnetic fields which do not depend on the charged particle attributes. One of 
the reasons for presenting the modeling sections in this paper are to give the reader an idea of the detail included in 
this modeling work. As a test case Hiatt and Wilbur’s (ref. 15) 9.2 cm diameter discharge chamber is modeled with 
the computational tool presented in this work. While we are not capable of presenting a comparison to Hiatt’s and 
Wilbur’s experimentally determined plasma ion production cost, it is the intent of the authors to be able to do this in 
the near future. Many different types of results for this 9.2 cm diameter discharge chamber are given in the Results 
section. Also primary electron number density results obtained with PRIMA are given so that the present results can 
be compared to them. 


II. Static Field Models 


A. Electromagnetic Fields 

The steady state magnetic field produced by the permanent magnets located on the discharge chamber walls can 
be modeled using a subset of Maxwell’s equations. The governing equation for a 2-D, axisymmetric magnetic field 
configuration is determined by looking at the magnetostatic part of Maxwell’s equations. This magnetic field can be 
determined with 


d_ 

dr 


1 


0(rTe)Y 


rji dr ) dz 


1 SA 0 
p dz 


8H cr dH cz 


dz 


dr 


(1) 


For the axisymmetric case of a 2-D model, the magnetic vector potential, Aq , only has a circumferential component. 
The axial and radial magnetic fields are related to the circumferential component of the magnetic vector potential as 


B = ! d(rAg) 


and 


B - = 


r dr 
dA e 


dz 


(2) 

( 3 ) 


Equation (1) is solved using following boundary conditions; 


Aq — > 0 as r — > co , 


Aq — > 0 as z — » -oo , 


( 4 ) 

( 5 ) 
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Aq — ^ 0 as z — ^ oo , 


(6) 


and 


9Aq 

dr 


= 0 at r = 0 . 


(7) 


To utilize these boundary conditions the static magnetic field computational domain must be made larger than the 
discharge chamber. This is done in this work. 

B. Electric Fields 

To determine the electrical potential caused by static electrical potentials, 4> applied to boundaries Laplace’s 
equation (ref. 16) is used 


V 2 ()) = 0. 


(8) 


The electric field produced by these static potentials is obtained from 


E = -V<|> . 


(9) 


The discharge chamber wall boundary conditions on the electric potentials are given by, 


* 1 * walls ~ 59 volts , 


(10) 


^cathode ~ 9 volts , 


(11) 


and 


§ screen 9 volts . 

III. Particle Models 


(12) 


A. Dynamic Fields 

For rapidly changing currents and rapidly changing charge densities a more sophisticated model of the electric 
and magnetic fields caused by the particles in the discharge chamber needs to be considered. This particular model 
will be called the dynamic electrical and magnetic field model. The dynamic electric and magnetic fields are 
modeled using the complete set of Maxwell’s equations. The integral form of Maxwell’s equations is written as, 
(ref. 16) 


Jf-d5 = JpdF, (13) 

s 
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(14) 


J5rf5=0, 

5 


o E ■ dl = - ^d t B dS , 


(15) 


and 


oH.dl = 1 8 t D dS + ^J-dS. 


(16) 


where E is the electric field, B is the magnetic field, H is the magnetic field strength, D is the electric 
displacement, p is the charge density, J is the current density, V is the volume, S is the surface area, / is the 
position vector and t is the time. Because of the principle of superposition for both the electric and magnetic fields 
the portion of the fields caused by static phenomena may be added to that caused by the dynamic phenomenon. The 
total electric and magnetic field at a given point in the discharge chamber is then obtained by adding the two effects. 
Utilizing this type of technique provides possibilities for speeding up the calculation. 


B. Particle Advance 

The governing equations for the particle motion are given by the classical Newton-Lorentz equations of motion 
(ref. 17). The charge particle equations of motion in the electric and magnetic fields are: 


nij 


dvj 

dt 


= <7/ E + V( x b\, 


(17) 


and 



(18) 


where w, is the mass of a particle, v ; is the particle velocity, qj is the electric charge of the particle and Xj is the 
particle position. The particle’s velocity and position are advanced along a trajectory from the point of its creation 
until the point it is lost by absorption to a boundary. In this work the electrons inside the discharge chamber are 
grouped in to two categories. They are: 1) primary electrons and 2) secondary electrons. Primary electrons are the 
electrons that are emitted from the cathode and they generally have higher kinetic energy than the secondary 
electrons. Secondary electrons are the electrons produced from the ionization collisions of electrons with the neutral 
atoms. Essentially both electrons are the same species but they are grouped separately in the analysis. The primary 
electrons are more important for the discharge chamber performance since they have higher kinetic energy and have 
more chances of producing ions than the secondary electrons. Modeling them as two groups is useful for 
presentation purposes. 

The boundary conditions for the particles are the particle emission, reflection and absorption. In a general ion 
engine discharge chamber the primary electrons are produced from a hollow cathode. The plasma near the cathode 
surface is quasi-neutral (ref. 18). At the exit of the hollow cathode orifice both primary electrons and ions are 
present. The ions are produced near the cathode tip mainly because of ionizing collisions of primary electrons with 
the neutrals. This ion production at the cathode tip facilitates a high electron emission current (ref. 19). In this model 
a plasma source boundary condition is used for the cathode emission. The rate of plasma, h p , produced is based on 
the emission current from the cathode, Ip , and the volume near the cathode tip, V c tip , chosen for plasma 
generation. It is given by 
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(19) 


hp 


Ie 

V(: tip 


where e is the absolute value of the charge of an electron. Both primary electrons and ions are produced at this rate 

at the cathode tip. The speed at which the primary electrons and ions are leaving the cathode tip is based on the 
discharge voltage, V d , applied and the cathode tip temperature respectively. The speed of the primary 
electrons, v pe | , leaving the cathode is determined using, (ref. 20) 


v 


pe 



(20) 


where m e is the mass of an electron and Vj is the discharge voltage in eV. Similarly the speed of the ions, |v, 0 „ 
that are leaving the cathode tip is calculated using, (ref. 20) 


3 W,. 


( 21 ) 


where m lon is the mass of an ion and T jon is the temperature of the ions at the cathode tip in eV. 

The electrons are reflected off of cathode biased walls and are absorbed by anode biased walls. The electrons are 
reflected specularly (ref. 20) back into the discharge chamber at cathode biased surfaces. The ions are absorbed at 
the discharge chamber walls and at grid locations. Generally the discharge chamber walls are set as anode biased 
walls and the screen grids at the downstream end of the discharge chamber act as cathode biased walls. All particles 
are specularly reflected back into the discharge chamber at the axis of symmetry boundary (ref. 1). 

C. Particle Collisions 

Particle collisions modeled for the electrons and ions are collisions with the neutral particles. Coulomb collisions 
such as electron-electron or electron-ion collisions are not specifically modeled in XOOPIC. In general, a PIC 
technique models the long range coulomb collisions and smoothes the short range coulomb collision effects 
(ref. 21). These collisions are modeled inherently in a PIC technique since the charge density at a given location is 
estimated based on the particle charges that are present within a grid cell volume. This contributes to the long range 
coulomb collisions, but the collisions of charged particles within a cell volume are underrepresented. Recombination 
of charge particles is not modeled in this work. 

The types of electron collisions with the neutral particles included in this model are elastic, excitation, and 
ionizing collisions. In an elastic collision, the electron loses little or none of its kinetic energy and gets scattered 
after the collision. In an elastic collision the electron can pass some of its energy of motion to the motion of the 
particle with which it is colliding even though none of the energy goes into increasing the internal energy of the 
particle. In an excitation collision, the electron loses some kinetic energy and gets scattered after the collision. The 
electron has to lose at least the excitation threshold energy. The excitation threshold energy is the amount of energy 
required to excite a neutral atom to the next energy level. In an ionizing collision the electron ionizes the neutral 
atoms and creates a secondary electron. The electron which caused the neutral to be ionized loses kinetic energy 
equal to or greater than the amount of energy required to ionize the neutral atom. This energy is known as the 
ionization threshold energy. Ion collisions with neutral particles are modeled as elastic or charge exchange collision 
types. In an elastic collision, the incident ion loses some or none of its energy and gets scattered. In a charge 
exchange collision, the incident ion excites a colliding neutral particle and after the collision the neutral particle 
becomes an ion and the incident ion becomes a neutral particle. 
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The occurrence of a particle collision is determined using a null collision technique (refs. 22 and 23). This 
technique uses a constant collision frequency to determine the maximum fraction of particles that undergo 
collisions. This method eliminates the look up of kinetic energy for each particle. The null collision probability 
(ref. 23) is given by, 


Pmdl = 1 - exp(-v c * At) , 


(22) 


where v c is the constant collision frequency obtained by looking for the maximum sum of all collisions over all 
possible energies. 


v c — max(n„ o j (8,„ c )|v ; „ c ) . 


(23) 


This is a computation convenience that is made to represent reality by adding a null collision frequency to the actual 
collision frequency at each energy level (see reference [23]). Here the neutral number density, n„ is given by 
(ref. 10) 


/?„ = [4m(l - t\ u )\/\e\v n \A g § n , (24) 

where A g is the grid area, m is the propellant mass flow rate, r\ u is the propellant utilization parameter, |v„ | is the 
velocity of the neural particles, and <j>„ is the grid transparency to neutral atoms. The total collision cross-section, 
a r (E,„ c ), is given by 


) — ^1 i l: inc ) + ^2 (S/'nc) + • • • + &N (fiinc') 


(25) 


where N is the total number of collision types considered. The collision cross section area for all collision types are 
based on the kinetic energy of the incident particle which is given by, 


&inc 


1 

^ Mine 



(26) 


where OT i; , c is the mass of the incident particle and |v,- nc | is the speed of the incident particle. Since uniform neutral 
number density is assumed the constant collision frequency becomes, 

v c = n„ * raax(G r (Sinc)\v in c\) ■ (27) 

The null collision probability is estimated using equation (22) with equation (27). The number fraction of particles 
that experience collisions is calculated based on this null collision probability, 

N coll = N total * Pmill ■ (28) 


Only this number of particles is tested for collisions and some of these particles may not undergo a collision. 


NASA/CR— 2005-213833 


7 




Figure 1. — A Flow chart of PIC-MCC model shows Figure 2. — Field variables are stored 

computing sequence of a plasma simulation problem as used in the Yee mesh. 

(from Vahedi and Surendra, 1995). 


IV. Solution of Mathematical Model 

In this section the solution of the mathematical models presented above is discussed. Of course numerical 
techniques are utilized to solve these equations. A general PIC-MCC (ref. 21) flow chart (see fig. 1) shows the 
computing sequence of how XOOPIC models a plasma simulation problem for each time step, At . In the following 
sections the details of how the field and particle calculations are advanced through time are given. Details on how 
the collisions fit into this process are also given. 

A. Electromagnetic Fields 

The field effects of any plasma device in XOOPIC are modeled on a two dimensional computational mesh. It 
handles orthogonal meshes of x-y coordinates and r-z coordinates. In our model we use cylindrical r-z coordinates 
since most of the ion engine discharge chamber configurations studied are cylindrically shaped. In XOOPIC the 
field variables, as shown in figure 2, are stored at the cell centers which are followed in the fashion as used on the 
Yee mesh (ref. 24). 

The electric and magnetic field problems are solved by utilizing two modes: static mode and dynamic mode. The 
static mode is solved once at the beginning of the program and utilized throughout the computation. The dynamic 
quantities are solved for each time step because they depend on the position of the particles. These two portions of 
the field solver can be handled separately because of the principle of superposition. 

The static magnetic fields are solved utilizing the computer code called MAXWELL-2D (ref. 25). MAXWELL- 
2D is a computer program developed by the Ansoft Corporation, Inc. A previous study by Menart (ref. 26) indicates 
that Maxwell-2D is an excellent software package for modeling two dimensional magnetic field problems. The 
details of the discharge chamber configuration, such as the location and orientation of the magnets and the magnet 
properties, are required inputs to Maxwell-2D. If the wall material interacts with the magnetic field their location 
and magnetic properties are also required. Maxwell-2D uses a finite element technique with an iterative or direct 
matrix solver to find the magnetic vector potential values. The magnetic field values are obtained using equations 
(2) and (3) which are then provided as an input to XOOPIC. 

The static electric fields are determined inside XOOPIC itself. XOOPIC has the option of choosing different 
types of Laplace equation solvers for any given problem. The Laplace solvers available in XOOPIC are: multigrid 
method, conjugate gradient algorithm, dynamic ADI algorithm, and generalized minimum residual algorithm. Any 
one of these algorithms can be used to solve Laplace’s equation. The potential values are obtained at the cell nodes. 
Then the corresponding E field values are calculated from equation (9). 

The dynamic portion of this field analysis is the most difficult. As can be seen in equations (15) and (16) the 
electric fields produced by time varying magnetic fields and the magnetic fields produced by time varying electric 
fields are included in the analysis. For the discharge chambers that are studied in this work neither one of these 
phenomena is very important. The currents produced in a discharge chamber are too small to produce any significant 
magnetic fields and charge densities do not change fast enough to cause electrical fields changes to produce 
noticeable magnetic fields. The important part of this dynamic analysis for ion engines is tracking the change in the 
electric fields caused by the changing particle locations. This is effectively done in XOOPIC. A set of magnetic field 
and electric field equations is obtained by using a time centered finite differencing on equations (15) and (16). 
Leapfrog scheme is adopted for the time integration of these field equations. The details of advancing 
electromagnetic field problems are discussed at length in reference 1 . 
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B. Particle advance 

XOOPIC uses the relativistic form of the equation of motion for charged particles in electric and magnetic fields. 
Time integration of these equations are handled using a Boris advance technique (ref. 17). The time centered finite 
difference form of the relativistic equation of motion is given as, 


■ 72 + 1/2 


— Uj 

At 


72 - 1/2 


- 72 + 1 /2 - 72-1 /2 
Qi - U; +Ui - 

— E n +— l - xB>‘ 

m i 2y i n 


(29) 


where Uj is the particle relativistic velocity which is equivalent to Uj = y/iy . The relativistic factor, y ,■ , is calculated 

I- |2 

, \Uj 

using the relation, yf = 1 + J — 2— ; where c is the speed of light. In the Boris advance the electric and magnetic forces 

c 1 

are separated in the integration by substituting the relations 


- n— 1/2 — qjE” At 

U; = Uj - - 


2m j 


(30) 


and 


Uj 


72 + 1/2 


= Uj 


qjE" At 
2 nij 


(31) 


into equation (29) giving a new time integration equation which eliminates the electric field forces from the finite 
difference equation. 


Uj - Uj 

At 


— ( m ; + +«/ ) *B n 


ry 72 

2y i m- 


(32) 


Equation (32) represents the rotation caused by magnetic field forces on the charged particles. The magnetic field 
rotation is handled by a Boris rotation which introduces another set of variables, 


Uj = Uj +Uj X tj 


(33) 


and 


Uj + = Uj 


+ Uj X Sj 


(34) 


with 


- _ qi B"At 
2y i'mj 


(35) 
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and 


S: 


(36) 


2 t t 


1 + tj 


The computational steps are: first compute the uf velocity by adding half the electrical force to the known 
velocity, at time level n- 1/2 using equation (30). Then rotate the uf velocity using equations (33) to 

(36) for a full time step to calculate z7, + velocity and add the remaining half electrical impulse, equation (31), to 
obtain the new velocity, m,-”' 1/2 at time level n + 1/2 . The particle positions are updated with new velocity. 


- n+ 1/2 . 

— n + 1 — n l\t 

X- = X, H 


Yi 


71 + 1/2 


(37) 


where y/ !+1 is evaluated at time level n + 1/2 . 

In an actual hollow cathode, the particles emitted from the cathode tip area, are from any arbitrary location 
within this area and with a random velocity. This random emission of primary electrons and ions from the cathode is 
handled using a Markov-Chain Monte Carlo scheme. The number of particles that are emitted for a given time step 
is determined with the source rate, cathode tip volume as given in equation (19). The speed of each primary electron 
and ion are calculated using equations (20) and (21) which utilizes the isotropic temperature values given as input. A 
Maxwellian velocity distribution is applied in each direction to get a velocity for each particle. For the primary 
electrons their temperature is equivalent to the discharge voltage of cathode. For the ions their temperature would be 
the temperature near the hollow cathode tip. 

C. Particle collisions 

The neutral particles are assumed distributed uniformly throughout the discharge chamber. The neutral 
temperature used in this model is 0.025 eV which is equivalent to 290 K. The neutral number density is calculated 
using equation (24). The maximum value of <^T( s inc)\^inc\ * s estimated by calculating the total collision cross 
section areas, Oj-(z inc ) , and the particle speed, \v- lnc \ , for a range of kinetic energies, z inc , from 0 to 1000 electron 

volts. From this range of particle kinetic energies the maximum value of the collision swept volumetric rate is 
calculated. Then utilizing the neutral number density and the maximum collision volume rate, the constant collision 
frequency, o c , is obtained using equation (23). The probability of null collision is calculated using a constant 
collision frequency and time step, At . Then the fraction of particles, N co u that undergoes particle collisions is 
estimated. The selection of which particles should have a collision is done through a random process. In this process 
effort is made not to duplicate the same particle for another type of collision in a given time step. Then the collision 
type for each colliding particle is tested using a random number, R , between 0 and 1 in the following manner, 

R < Ui(e,„ c )/o c (Collision Type 1) 

Ui(e inc)/ V C <R< (ui(e,-„ c ) + o 2 (s inc))/ o c (Collision Type 2) 


N 

7o j(£j nc )/ u c < 7? (Null collision) (38) 

7=1 


NASA/CR— 2005-213833 


10 



This collision process is performed for the primary electrons, the secondary electrons and the ions in the same way. 
In the case of electron particles three collision types: elastic, excitation and ionization are modeled. For ion particles 
two collisions types, elastic and charge exchange, are modeled. The collision model used in XOOPIC to handle 
these particle collisions are from the work of Vahedi and Surendra (ref. 23). More details about collisions can be 
obtained from their paper (ref. 23). Also XOOPIC uses collision cross-section areas for the elastic and excitation 
collisions from Surendra et al. (ref. 27) and the ionization collision cross section areas from Rapp and Golden 
(ref. 28) for the argon gas. 

D. Stability Criteria 

XOOPIC uses an explicit finite difference time domain (FDTD) algorithm for solving the electromagnetic field 
equations and an explicit leap frog time mover for the particle advancing. These explicit schemes have limitations 
on selecting the grid spacing sizes, A z and A r , and the time step size, At . The FDTD algorithm for the 
electromagnetic solver needs to meet the Courant condition (ref. 17) to control the nonphysical growth of numerical 
errors. The Courant condition is 



1 

A r 2 


< 1 , 


(39) 


where c is the speed of light. The explicit particle advancing scheme (ref. 17 and 21) is limited by the following 
stability constraint on the time step size, 


?n p e At A 2 , 


(40) 


and for accuracy. 


tn pe At « 0.2 , 


(41) 


where m pe is the plasma frequency and it is given by (ref. 29) 


I |2 

n e \e\ 

V £o m e 


(42) 


The constraint for grid spacing is given as, (ref. 21) 


^ D 


1 


1 


A A r~ 


>0.3, 


where Xj) is the Debye length and it is given by, (ref. 29) 


(43) 
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Here n e is the electron number density and T e is the electron temperature. The constraint for grid spacing is 
required to control the plasma heating (ref. 17) which happens with larger grid sizes. The particle weighting is also a 
factor to maintain the stability of the numerical plasma simulations. The numerical noise is proportional to ^jl/N cp , 
where N cp is the number of macro particles in the simulation. Weighting each particle with fewer physical particles 
will reduce the numerical noise. In general too few particles can result in numerical heating. 

V. Test Model 

In this work an experimental discharge chamber design 
studied by Hiatt and Wilbur (ref. 15) is used as a test model. 

The reason for choosing this model is because previous 
modeling studies have used this chamber design for 
comparisons. Hiatt and Wilbur used a 9.2 cm diameter steel 
walled discharge chamber to perform experimental 
investigations. The chamber length used on this model is 
11.8 cm. Figure 3 shows the schematic view of this 
experimental model. It is a straight cylindrical shaped 
discharge chamber. It has two samarium cobalt magnets 
which are placed inside the discharge chamber. Magnet 1 is 
placed on the back wall and Magnet 2 is placed on the side 
wall at a distance of 3.7 cm from the screen grid. The 
magnetic flux density of the magnets is maintained at 2700 
Gauss. The propellant used in this model is argon. An anode 
plate is placed axially at the bottom surface of Magnet 2. A circular, filament type cathode of 1 cm diameter is used 
in this experimental setup. This cathode is made from two coiled tungsten wires of 0.25 mm diameter which are 
wound together to form a spiral hoop which has a diameter of 1 cm. This hoop shape keeps the cathode off the 
center line. This cathode is located 1.5 cm upstream of the screen grid surface. The electric potentials at the wall 
boundaries, anode surface, cathode surface and screen grid surface need to be specified before performing 
simulations in XOOPIC. From Hiatt and Wilbur’s paper it is found that screen grid was maintained at 750 volts and 
the discharge voltage was 50 volts. Generally in ion engine discharge chambers, the cathode and the screen grid 
surfaces are maintained at the same potential (ref. 30) and the anode walls are kept at a higher potential to get the 
cathode discharge voltage. This procedure is followed in this work. The cathode, the discharge chamber walls, and 
the screen grid are at zero potential and the anode plate is maintained at 50 volts. The neutral number density is 
found from the relation given by equation (24). This equation gives a neutral number density of 3.3x10 (ref. 19) 
#/m 3 for the mass flow used in this paper. 

Hiatt and Wilbur performed tests for a range of argon mass flow rates from 100 to 700 mA eq. They measured 
the discharge current. Ip , , the ion beam current, Ip , and the total ion production current. Ip for each mass flow 
rate. These measured currents were used to compute the plasma ion production cost, e p, and the extracted ion 
fraction, fg utilizing the relations below (ref. 15): 



Figure 3. — Hiatt and Wilbur’s experimental model. 


Ib = Ib! Ip, 


(45) 


and 


Ep = ( Id ~ Ip) v d! Ip ■ 


(46) 
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Hiatt and Wilbur present e p results in a plot against the neutral density parameter, m( 1 - r\ u ) , for each mass flow 
rate. In this work the case of maximum mass flow rate, 648 mA eq., and the neutral density parameter, iii(\ - r \ u ) , of 
610 mA eq is taken. For this test case, the plasma ion production cost, £p is found to be 70 eV from their plot. Also 
the extracted ion fraction, fg is found to be 0.23 from their paper. Also needed are the cathode emission current, 
/ p . and the cathode tip volume to calculate the plasma source rate. From Brophy’s model, the cathode emission 
current is given by, (ref. 10) 


Ie =1 D~Ip ■ 


(47) 


Substituting sp and fg values into equations (45) to (47) gives the emission current. The emission current is found 
to be 232 mA. The cathode tip volume is determined by assuming the tip as a circular cylindrical ring. It is given by, 

Vc_tip = n(r?-r£ lt )h (48) 


where r- is the inner radius of the cathode hoop, r out is the outer radius of the cathode hoop, and h is the length of 

the cathode coil. Since the cathode details from Hiatt and Wilbur’s paper were not clear we assume a 1 mm circular 
coil thickness and a 1 mm length for the coil. This assumption gives an inner hoop radius of the cathode wire of 
5 mm, an outer hoop radius of 6 mm and wire width of 1 mm. Substituting these values into equation (48) gives the 
cathode tip volume. Using equation (19) the plasma source rate can be obtained. The plasma source rate used in this 
work is 3.6x10 (ref. 25) #/(m 3 -s). This parameter is used to get the number of primary electrons and ions emitted 
from the cathode tip at a given time step. 

The numerical parameters chosen to perform this simulation in XOOPIC are: the grid spacing, A z and A r , the 
time step, At, and the particle weighting. These parameters are given in equations (39), (40), (41), and (43). The 
grid spacing parameters are based on the Debye length and the time step size is based on the plasma frequency. X p 
and m pe need to be estimated for this discharge chamber model. For the plasma in the discharge chamber of an ion 

engine a high value of the electron number density (ref. 31) is lxl0 18 #/m 3 . This number density value with a range 
of electron temperatures from 10 to 50 eV, is used in equation (44) to obtain a set of Debye length values. The 
axisymmetric chamber dimensions of 11.8 cm length and 4.6 cm radius are used to get the grid spacing values. The 
grid spacing values and the debye lengths are checked to select a grid spacing that meets the debye length limit. 
Using 1000 grid points in the axial direction and 500 grid points in the radial direction satisfies the debye length 
limit for the range of electron temperatures mentioned above. Similarly the plasma frequency is calculated using 
equation (42) by utilizing the same high electron number density value. It is found to be 5.63xlO ln s _1 . The Courant 
condition in equation (39) puts a smaller limit on the time step value. For this reason a time step value of 2xl0~ 12 s 
was used. This value worked well. We can select an even finer grid spacing and a smaller time step to improve the 
accuracy of these simulations. However, at present time the grid spacing is limited by the availability of computer 
memory. A smaller time step is not used to keep the computational time from becoming undoable. The particle 
weighting also plays a role in the accuracy. In this work a particle weighting of lxlO 6 real particles per 
computational particle is used to maintain the numerical noise small while not making the computations excessively 
long. The numerical parameters considered used in the computation are: 


A z = 1.18x10 4 m , 

(49) 

Ar = 0.92xl0 -4 m , 

(50) 

At = 2xl0~ 12 ^ 

(51) 

lxlO 6 real particles/computational particles . 

(52) 
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VI. RESULTS 

In this section the results from Maxwell-2D and XOOPIC for the test model are presented. From Maxwell-2D 
the magnetic field results are obtained. These results provide information about the magnetic field values at different 
locations inside the discharge chamber. These values are essentially constant for all times in the discharge chamber 
and thus they are only shown once. 

From XOOPIC numerous results are obtained. The results are presented in two categories: particle and electric 
field. For the most part these results are given at 2 and 300 ns, except for the electrical potential and the total number 
of particles in the discharge chamber. The electrical potential is shown at time zero so that the electrical potentials 
caused by the applied boundary conditions can be seen. Any change in these values is caused by the particles. The 
total number of particles is shown as a function of time in the discharge chamber. Results from 2 and 300 ns are 
presented next to each other so that comparisons can be made. The last result presented in this paper is that from 
PRIMA. This result gives the primary electron number density in the discharge chamber at steady state conditions. 


A. Initial Results 

The magnetic vector potential contours are given in figure 4. This plot shows how the magnetic field lines run 
between the two magnets. A 100 gauss-cm contour line forms near the cathode and ends at the center of the anode 
plate placed below Magnet-2. This contour line will tend to direct the motion of the primary electrons. Figures 5 and 
6 show the contour plots of radial and axial magnetic fields. Strong negative radial field lines are observed below the 
Magnet-2 surface. The magnetic field strength is decreasing radially as we move downward from the second magnet 
axial surface. Also looking from figure 5 and figure 6 we can see that the cusp region at the Magnet-2 surface has 
stronger radial field components and weaker the axial field components. This indicates that the radial moving 
electron particles will have an easy path of moving into the cusp region and axially moving electron particles will be 
confined by the strong radial magnetic field. 

The electric potential at time zero is shown in figure 7. This result is obtained from XOOPIC by solving the 
Laplace equation of electric potential. This plot shows zero potential value at the discharge chamber walls and at the 
screen grid areas. Also the potential value at the anode plate surface is 50 volts. This 50 volts value is diffusing to 
zero value around the anode wall surface. 



Zinin 


Figure 4. — Magnetic vector potential contour lines Figure 5. — Radial magnetic field contour lines in gauss, 

in gauss-cm. 




Zinin 

Figure 7. — Electric potential contour lines in volts. 
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B. Particle Results 

In this section the results of primary electrons, secondary 
electrons and ions inside the discharge chamber are presented. 
A plot of total number of primary electrons, secondary 
electrons, and ions present inside the discharge chamber is 
shown in figure 8. Both the vertical and horizontal axes in this 
plot have a logarithmic scale. The total number of particles 
includes all the physical particles present inside the discharge 
chamber. This 

value is obtained by multiplying the number of macro particles 
with its weighting. This plot indicates that all particles are 
steadily increasing with time. Primary electrons and argon ions 
are increasing at the same rate until 40 ns. After 40 ns, the ions 
are increasing at a slightly higher rate than the primary 
electrons. This difference in rate of increase for ions and 
primary electrons is due to the increase in number of ionizing 
collisions of primary electrons with the neutrals. The first 
secondary electron is formed at 0.8 ns. By looking at figure 8, 
we can say that steady state operation of this discharge 
chamber is not yet achieved and the simulation still needs to be 
run for a long time. It took 5 days of computational time on a 
single processor to compute this simulation up to a time of 
300 ns. This is the limit of computational time available for 
this work. In the future the computer code will be parallelized 
and steady state will be reached. 

The primary electron number density results at 2 and at 
300 ns are presented in figure 9 and figure 10 respectively. By 
Looking at figure 9 one can see that few primary electrons are 
present inside the discharge chamber at 2 ns. The ones that are 
present at 2 ns are located near the cathode. From figure 10, it 
can be seen that the region of primary electron movement is 
along the magnetic field lines that run from the cathode to the 
anode located above the side wall magnet (see fig. 4 and 
fig. 10). While it is difficult to see in figure 10 some primary 
electrons are making their way towards the cusp region of the 
back wall magnet. For the most part the primary electrons are 
spread near the cathode radially for a height of 1.2 cm. The 
maximum number density value of 3.3xl0 ls #/nr’ is observed 
above the cathode and the number density value near the 
screen grid is about 10 I4 #/m’. 

The secondary electron number density results are also 
interesting to see. These electrons are created due to the 
ionization collisions of primary electrons with the neutrals. 
Figure 11 shows the secondary electron number density 
contours at 300 ns. The results of 2 ns case are not shown here 
because the number of secondary electrons produced at 2 ns is 
very small, around 2xl0 7 . At 300 ns, a significant number of 
secondary electrons are present inside the discharge chamber. 
A high value of secondary electron number density, lxl0 ls 
#/m 3 is observed near the cathode location. The secondary 
electron number density decreases rapidly at locations away 
from cathode. Only a few secondary electrons are present in 
the cusp region at 300 ns. 

The argon ions are produced from the cathode tip and from 
ionizing collisions of electrons with the argon neutrals. The 
argon ion number density values inside the discharge chamber 
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Figure 8. — Total number of particles vs. time. 
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Figure 9. — Primary electron number density 
contour lines in #/m 3 at 2 ns. 
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Figure 10. — Primary electron number density 
contour lines in #/m 3 at 300 ns. 




Figure 11. — Secondary electron number density 
contour lines in #/m 3 at 300 ns. 
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Figure 12. — Argon ion number density 
contour lines in #/m 3 at 2 ns. 
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Figure 13. — Argon ion number density 
contour lines in #/m 3 at 300 ns. 


at 2 ns are given in figure 12. The maximum number 
density value is observed at the cathode. The speed of ions 
is two to three orders of magnitude less than the electron 
speed because ions are much heavier than electrons. This is 
the reason the ions have not moved as far as the primary 
electrons at 2 ns (see fig. 9). Figure 13 shows the argon ion 
number density values at 300 ns. This plot shows the ions 
moving away from the cathode region. The interesting 
result from this density plot is that the ions are moving 
radially much further into the stronger axial magnetic field 
region. This stronger axial field significantly confines the 
electrons, but the ions cross these lines and move further 
into the strong regions. This movement of ions into the 
stronger magnetic field region is not surprising because the 
Larmor radii (ref. 16) for the argon ions are much higher 
than the Larmor radii for the electrons. A high value of the 
argon ion number density of 4.5xl0 18 #/m 3 is observed at 
the cathode. 

C. Field Results 

In this section the results of the electric potential are 
presented. The electric potential results at 2 and at 300 ns 
are shown in figures 14 and 15. The electric potential 
values observed near the cathode area are high in both 
cases. From figure 14, it can be seen that the region of high 
potential is limited to the cathode tip area because all argon 
ions are confined in this small region. There is a net 
positive charge density in this region because the primary 
electrons move out quicker than the ions. Results at 300 ns 
shown in figure 15 indicate a wider area of high electric 
potential. This is due to the movement of ions away from 
the cathode tip into the interior of ion engine discharge 
chamber while the primary electrons continue to move 
away even further. 


0 20 50 100 
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Figure 14. — Electric potential contour lines 
in volts at 2 ns. 
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Figure 15. — Electric potential contour lines 
in volts at 300 ns. 



D. Comparisons 

To get a feel for how the computational results from XOOPIC are comparing to results obtained from PRIMA 
figure 16 is given. The results in figure 16 show the relative primary electron number densities as obtained from 
PRIMA for the Hiatt and Wilbur’s model. PRIMA is a computational tool that is available to model only the 
primary electron motions inside the discharge chamber. It does not include the other physics which are modeled in 
XOOPIC such as electric field effects, modeling of secondary electron, modeling of ions, etc. PRIMA gives two 
primary outputs as results: primary electron confinement length and primary electron relative number density. The 
results obtained from PRIMA for the Hiatt and Wilbur’s model have good comparisons (refs. 8 and 32) with Hiatt 
and Wilbur’s results for the primary electron utilization factor. 
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The relative primary electron number density, n , is determined using, (ref. 6) 


* 



n 


pe 


n 


pe, max 


(53) 


where n is the primary electron density at some location and n pe max is the maximum primary electron number 

density. The values of relative primary electron number density vary between 0 and 1. Generally the value of 1 is 
observed at the cathode tip. The maximum value of primary electron number density observed at 300 ns is 
3.3xl0 l8 #/m 3 from XOOPIC. This value can be used to convert PRIMA’s relative primary electron number density 
results into an actual primary electron number density distributions. This converted result is shown in figure 16. 

Comparing figure 16 with figure 10 we can see that the number density contours at PRIMA and XOOPIC follow 
the same basic trend toward the side wall magnet, but different towards the back wall magnet. More than likely this 
is due to the fact that XOOPIC results have not made it to steady state. In future work XOOPIC will be taken out to 
steady state and this comparison will be made again. At this time all that will be said is that it appears as if results 
from XOOPIC are beginning to look like those from PRIMA. 

At this time XOOPIC is not capable of giving primary electron confinement length so a comparison for this 
output can not be made. In the future XOOPIC will be modified to produce this quantity. 

VII. Future Work 

The results presented in this work are preliminary in the computational modeling of an ion engine discharge 
chamber. In this work it is demonstrated that using a computational tool called XOOPIC it is possible to model the 
primary electrons, secondary electrons and ions inside the discharge chamber including the electric and magnetic 
field effects. The results are studied on an experimental design of discharge chamber of an ion engine. In the future 
it is hoped to get the steady state results for this model. Currently XOOPIC is not running in a parallelized mode. 
Parallelization has to be utilized in modeling this discharge chamber design to reach steady state results. 

In addition, the Debye length limit poses a big problem in resolving the computational mesh, thus non uniform 
mesh model needs to be explored. 


VIII. Conclusion 

In this work it is demonstrated that primary electrons, secondary electrons, and ions can be simultaneously 
tracked throughout a discharge chamber while the detailed electrical fields are being updated to include the effects 
of the charged particles. In this work a number of results are presented for a 9.2 cm diameter discharge chamber that 
is 11.8 cm long. In these it is seen that the ions move into the stronger magnetic field regions more than the 
electrons. The primary electrons tend to move out away from the cathode much quicker than the ions or secondary 
electrons. The primary electrons tend to ran along magnetic field lines that run from the cathode to the anode. The 
secondary electrons tend to stay in regions where the ions are located. The electrical fields are greatly affected by 
the charges present at the cathode. In general the cathode region is the location of highest particle number density 
for all particles computed. While not shown in this paper a number of other results can be obtained from XOOPIC 
like the electric field, the average kinetic energy and the particle flux to the walls and to the anode surfaces. The 
computer program XOOPIC is capable of producing a great deal of information on the discharge present in an ion 
engine. 
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